function out = getQMLse(paramsValues,setupFilter,outFilter)

numParams  = length(paramsValues);
stepOption = 2;
[hessian,scoreQML,scoreSMM,dQ,dQMLpart,dSMMpart,errorMes] = getScoreHessian(paramsValues,setupFilter,outFilter,stepOption);
if setupFilter.lambdaSMM  == 0
    score    = -scoreQML;
else
    score    = -scoreQML + scoreSMM;
end
% For debugging
% sum(scoreQML/(setupFilter.T-setupFilter.numInitial),1)
% dQMLpart'
% sum(scoreSMM/(setupFilter.T-setupFilter.numInitial),1)
% dSMMpart'
% sum((scoreQML+scoreSMM)/(setupFilter.T-setupFilter.numInitial),1)
% dQ'

% The inverse of the Hessian
invHessian       = hessian\eye(numParams);
if errorMes ~= 0
    out  = NaN;
    return
end

% The variance of the score function using the Newey West estimator
qLag_index = [0 1 3 5 10 15 20];
out = struct();
for i=1:size(qLag_index,2)
    varScore                = NeweyWest(score(setupFilter.numInitial+1:end,:),qLag_index(1,i));
    AVar                    = 1/(setupFilter.T-setupFilter.numInitial)*invHessian*varScore*invHessian;
    name                    = ['qLag',num2str(qLag_index(1,i))];
    out.(name).paramsSE     = cell2struct(num2cell(sqrt(diag(AVar))),setupFilter.selectParams);
    out.(name).paramsTstats = cell2struct(num2cell(paramsValues./sqrt(diag(AVar))),setupFilter.selectParams);
end


end